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Abstract 

We present computational results showing eccentricity growth in the inner portions 
of a protoplanetary disc. We attribute this to the evolving surface density of the disc. 
The planet creates a gap, which adjusts the balance between the 3:1 (eccentricity 
exciting) and 2:1 (eccentricity damping) resonances. The eccentricity of the inner 
disc can rise as high as 0.3, which is sufficient to cause it to be accreted onto the 
star. This offers an alternative mechanism for producing the large holes observed in 
the discs of CoKu Tau/4, GM Aur and DM Tau. 
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1 Introduction 



Nearly two hundred extra-solar planets have been discovered over the past 
decade via the radial velocity (RV) method. Unfortunately, the RV method 
can only detect planets in tight orbits around older stars, making it unsuitable 
for detecting planets in the early stages of their formation. This deprives us 
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of key information about the processes of planetary formation, meaning that 
many of the most fundamental questions about the formation process remain 
unanswered. 



Recent observations by the Spitzer telescope have sho wn that the accretion 



discs surrounding s ome young stars contain large holes (ID'Alessio et al.l . 12005 



Calvet et all 120051 ). A variety of mechanisms for forming these holes has been 



proposed, but the most likely involve a planet 'damming up' the outer disc 
while the inner disc dissipates. The holes have sharp ed ges, and only a mas- 
sive planet provides a satisfactory explanation of this. iQuillen et al. (l2004t ) 
invok ed viscous diss i pation of the disc to explain the hole in CoKu Tau/4, 
while IVarniere et al.l (120061 ) suggested that waves launched from Lindblad res- 
onances might be responsible. In this paper, we present numerical results which 
suggest another possibility, that of eccentricity excitation. 



2 Numerical Setup 



We use the Fargo code of Masset ( 2000al lbl) to perform our calculations. 



Fargo is a simple 2D polar mesh code dedicated to disc planet interactions. 
It is based upon a standard, ZEUS-like hydrodynamic solver, but owes its 
name to the Fargo algorithm upon which the azimuthal advection is based. 
This algorithm avoids the restrictive timestep typically imposed by the rapidly 
rotating inner regions of the disc, by permitting each annulus of cells to ro- 
tate at its local Keplerian velocity and stitching the results together again 
at the end of the timestep. The use of the Fargo algorithm typically lifts 
the timestep by an order of magnitude, and therefore speeds up the calcula- 
tion accordingly. Farg o includes an artifical visc osity like that described by 



equations 33 and 34 of IStone and Norman! (119921 ). but with C2 = 2 (thereby 



spreading shocks over two grid zones). The mesh centre lies at the central 
star, so indirect terms coming from the planets and the disc are included in 
the potential calculation. We make use of an open inner boundary. If v r < 
in the innermost active cells, it is copied to the ghost cells. However, if v r > 
in the innermost active cells, then the radial velocity in the ghost cells is set 
to zero. 

We use units normalised such that G = + M p = 1, while the planet's 
initial orbital radius is set at r = 1. An a-type viscosity is assumed, with 
a = 10~ 4 for the physical viscosity. We assume a constant aspect ratio disc, 
with h/r = 0.05, making the viscous timescale 
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approximately 4 x 10 5 orbits at r = 1. We set the mass ratio q = M p /M^ to be 
10~ 3 , approximately equal to the Jupiter-Sun value. The initial disc surface 
density is X oc r _1 , a reaso nable 'in t ermed iate' slope which makes the disc- 
only solution stationary (see iPringld . 119811 ) . The density is set to 1 x 10~ 5 at 
the planet's radius in our units. This means that, although the planet is free 
to migrate, there is insufficient mass present in the disc to cause significant 
orbital evolution of the planet. Our introduction of a planet into a smooth 
disc implicitly assumes that the planet has gone from the gap opening criteria 
to a Jupiter mass ve ry rapidly. This is reasonable in bo t h the core accretion 



( [Pollack et al.l . Il996l ) and gravitational instability (jBossl . 120031 ) models. 



Our computational domain extends between r = 0.07 and r = 4. It is covered 
by 128 grid cells in the radial direction (logarithmically spaced) and 384 az- 
imuthally (uniformly spaced). The effect of the planet's gravitational force is 
smoothed over 0.6 of the disc thickness. 



2. 1 Test without a planet 



As an initial test, we ran the code without a planet present, but with all other 
parameters as described above. Figured] shows the evolution of inner disc mass 
and eccentricity with time, while figure [2] gives three snapshots of the surface 
density. As expected, the disc is almost stationary. On examining the viscous 
time scale as a func tion of radius, we found it to follow the standard prediction 
(e.g. |Pringld . [l98lh . 



2.2 Eccentricity in the Outer Disc 



Kiev and Dirksenl (l2006h recently presented results showing a planet inducing 
growth of eccentricity in the outer disc. We have run two comparisons to 
th eir code, with a q = 10 -3 and q = 4 x 10~ 3 planet. According to the work 
of iKley and Dirksenl . the Jupter mass planet should not induce significant 
eccentricity in the outer disc, whereas the 4 Mj planet should. 

Figure[3]plots the radial kinetic energy of t he disc for the two planetary masses. 
It should be compared with figure 8a of IKley and Dirksenl . We see that the 
4 Mj planet has excited the eccentricity, whereas the 1 Mj planet has not. In 
figure HI we plot the radial eccentricity profil e of the discs, afte r 2500 orbits. 
We see that our results are similar to those of IKley and Dirksenl . In a separate 
test with a grid e xtending to r = 2 .5, we found our peak eccentricity to be 
suppressed, as did lKley and Dirksen . 



Although the value of eccentricity reached in our calculation is similar to that 
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Fig. 1. Evolution of the mass and mass averaged eccentricity of the inner disc 
(0.07 < r < 1) without a planet. The eccentricity is calculated by assuming that 
each grid cell is a free particle 




Fig. 2. Evolution of the disc surface density in the inner disc when no planet is 
present. The units of surface density are those described in the text 



observed by lKley and Dirksenl (120061 ). the eccentric mode grew more rapidly in 
our calculation. The radial kinetic energy in our test peaks after around 1000 
orbits in our test, but after 1500 orbits in the the work of iKlev and Dirksen 



(2006). We believe that this is due to the effect of the disc mass itself. iLubow 



1991al ) predicted that the eccentricity should grow on a timescale proportional 
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Fig. 3. Evolution of r adial kinetic energy of the disc, as a comparison with 
Kiev and Dirksenl (|200d ). Two different planet masses are plotted 
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Fig. 4. Radial eccentr icity profiles after 2500 orbits, as a comparison with 
Kiev and Dirksen tod ) for two different planet masses 
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where the surface density, E, is evaluated at the location of the exciting 
resonance. In a linear theory, the density should be proportional to the to- 
tal disc mass, and hence the disc mass itself should be irrelevant. However, 
we are not in the linear regime, and equation [2] takes no account of the 
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Fig. 5. Evolution of the mass and mass averaged eccentricity of the inner disc 
(0.07 < r < 1) with a q = 10~ 3 planet. The eccentricity is calculated by assuming 
that each grid cell is a free particle 

damping resonanc e s, or how the eccentric mode propagates through the disc. 



Kiev and Dirksenl ( 120061 ) imply that their disc mass was actually on the or- 
der of the stellar mass, which is unrealistically high. In further tests, we have 
found that increasing the disc mass can slow, or even suppress, the growth of 
eccentricity in the disc. 



3 Results 



Having tested FARGO, we now proceed to our main calculation: that of eccen- 
tricity in the inner disc. 

In figure [5] we plot the evolution of the mass and mass averaged eccentricity 
of the inner disc (between r = 0.07 and r = 1). The eccentricity is quite 
low for around 500 orbits, when it suddenly starts to climb rapidly. As its 
eccentricity grows, the inner disc is depleted (although there is a slight lag 
behind the eccentricity growth). We attribute this depletion to material falling 
onto the central star (i.e. falling off the computational grid). The eccentricity 
peaks after around 1500 orbits, and then decays. The outer disc stays in a low 
eccentricity state, with e < 0.02. The depletion of the inner disc is far faster 
than the viscous timescale computed above. The eccentric inner disc undergoes 
retrograde precession on a timescale of roughly 80 orbits. We present a sample 
snapshot of the density field in figure [61 The eccentricity of the inner disc is 
plain in this image. 
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Fig. 6. Snapshot of the density after 1500 orbits with a q = 10 3 planet. The density 
is plotted relative to the initial value, on a linear scale extending [0, 1] 

Th e most l ikely e x planatio n for this growth lies in the mechanism explored 
by iLubowl (1991aj). iLubowl describes how a circumstellar disc in a binary star 
system can become eccentric due to interactions with eccentric Lindblad res- 
onances. In particular, the eccentric Lindblad resonance located at r = 0.5 
(the 3:1 mean motion resonance) does not concide with an eccentricity damp- 
ing corotation resonance. Instead, any eccentricity excited must be damped 
through the 2:1 corotation resonance (located close to r = 0.6). These reso- 
nances are finely balanced, and their strengths determined by the local surface 
density. In figure [7] we show how the surface density evolves in our calculation. 
Between the start of the calculation and 500 orbits, we see that the density at 
the eccentricity exciting 3:1 resonance (r = 0.5) increases substantially, while 
the density at the eccentricity damping 2:1 resonance (r = 0.6) increases by a 
smaller fraction. This relative increase in excitation enables the 3:1 resonance 
to 'win' and excite the disc eccentricity. Material is then lost into the hole in the 
centre of the grid. As this happens, the density at the 3:1 resonance becomes 
comparable to that at the 2:1, and damping increases once more. Plotting 
eccentricity as a function of radius, we find that the eccentricity peaks close 
to r = 0.3, not at the 3:1 resonance itself. This lends support to the finding 
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Fig. 7. Evolution of the disc surface density in the inner disc when a q = 10 -3 planet 
is present. The units of surface density are those described in the text. The density 
increase at the inner edge of the disc for the 1500 orbit curve is confined to the two 
innermost grid cells 



of iGoodchild and Ogilvid ( 120061 ) that eccentricity could be suppressed at the 



location of the exciting resonance. 

In figure [HI we present a set of four runs, each showing growing eccentricity 
depleting the inner disc. These vary the initial surface density profile (using 
S oc r _1 and constant S) and the inner boundary condition. Two use the 
open inner boundary, described above, while two make use of Fargo' s 'non- 
reflecting' boundary. This attempts to match the guard cells smoothly onto 
waves propagating on the grid, so that no waves are reflected at the boundary. 
We see that in all cases, the eccentricity grew, and depleted the disc. Changing 
the boundary had little effect, while changing the intial surface density caused 
the eccentricity to appear earlier. A constant surface density is not a steady 
solution of the viscous disc evolution equation, so viscous evolution will have 
been occuring simultaneously with the gap opening and eccentricity excitation. 

Figure [9] shows the evolution of a disc identical to that in figure [5] except 
that it had a = 10~ 3 (that is, ten times greater than the assumed viscosity 
for figure [5]). We see that the disc undergoes more viscous evolution (as one 
would expect), and that the onset of eccentricity growth is delayed. 



Kiev and Dirksenl (120061 ) did not see eccentricity growth in their inner discs 
because their computational domain was too small. Their grid only extended 
to r = 0.25, which limits the eccentricity which could develop before material 
would stray from the grid. Furthermore, their boundary conditions strongly 
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Fig. 8. Evolution of eccentricity and inner disc mass in four runs with a q = 10 
planet. Top row: Open inner boundary; Bottom row: non-reflecting inner boundary 
Left column: £ oc r _1 initially; Right column: £ initially constant. Hence, the figure 
in the top left is identical to figured] 




Fig. 9. Evolution of the mass and mass averaged eccentricity of the inner disc 
(0.07 < r < 1) with a q = 10~ 3 planet and a increased to 10 -3 . All other pa- 
rameters were identical to figure [5] 
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suppressed eccentricity growth in the inner disc. iKlev and Dirksenl use d 'wave 



killing' boundary conditions (also used in Ide Val-Borro et all 120061 ) which 



push the solution within a 'wave killing' region back towards the standard 
Keplerian disc solution. Their boundary region extended to r = 0.5, the lo- 
cation of the 3:1 exciting resonance. Hence, any eccentricity excited at this 
resonance would have been subject to stro ng numerical damping. E ccentricity 



excitation of a protoplanet was studied by IPapaloizou et al 



(1200 If ), but they 



did not consider the effect on an inner disc. In a separate test, we found that 
a grid which only extended to r = 0.2 did not show such strong eccentricity 
growth (although computational speed was greatly increased). Instead, the 
material in the disc drained steadily, while remaining in near circular orbits. 

Explaining the unequal clearing of the resonance s is not completely stra ight- 
forward. In conventional theory (e.g. the review of lWard and Hahnl . |2000| ) . the 
disc interacts with the planet at its Lindblad resonances (LR). Most workers 
assume local damping, where the angular momentum in the excited wave is 
immediately transferred to the disc. In the inner disc, the innermost LR is the 
m = 2 which is colocated with the 2:1 resonance. This would appear to offer 
an easy explanation: the m = 2 LR clears material away, weakening eccen- 
tricity damping. There is no similar clearing from the 3:1 exciting resonance, 
since all other LR lie closer t o the planet. Howeyer, in common with many 



other workers in this area (see Ide Val-Borro et all 120061 ) , we see the planet 's 



wake propagating far into the inner disc, indicating that damping is not com- 
pletely local. Some of the wave may be damped immediately upon launch, 
but a noticeable fraction remains to propagate through the disc. It is possi- 
ble that a sufficient portion of the wave is damped at the resonance to clear 
the 2:1 damping resonance, but we cannot state this with certainty from our 
present calculations. We note that the imposition of a locally isothermal equa- 
tion of state p revents waves from steepening into shocks as they propagate (cf 
Rafikovl - liooi ). 



4 Discussion 



Disc eccentricity in an inner disc has been invoked to e xplain the 'super- 
hump' phenomenon in binary stars (e.g. iTruss et all 120011) . Indeed, this was 



the original motivation for the work of Lubow ( 1991a ,bf) Scaling typical re- 
sults according to equation [2], we find that the results presented above have a 
surprisingly rapid growth rate for the eccentricity. However, there are reasons 
why a straightforward scaling may not be appropriate. Firstly, there is the 
issue (mentioned above) of the disc mass. Our results are certainly affected 
by the disc mass, while equation [2] (based on a linear calculation) predicts 
that this should drop out of the growth rate. However, a Jupiter mass planet 
opens a gap in the disc on a timescale of 100 orbits or less, thereby putting the 
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system into the non-linear regime. Secondly, most work on binary stars uses 
the smoothed particle hydrodynamics (SPH) scheme to model the behaviour 
of the disc. SPH is known to possess a high numerical dissipation, and this 
could well affect the growth rate of the eccentric mode. The differences be- 
tween grid-based and SPH- based calculations o f plan ets in discs is vividly 
demonstrated by figure 10 of Ide Val-Borro et al.l (120061 ). Being performed over 
a decade late r, these SPH codes used nearly two orders of magnitude more 
particles than lLubowl (albeit over a larger disc), and yet still obviously have far 
higher numerical dissipation than the grid based codes. In figure EE we have 
shown how increasing the viscosity delays the onset of eccentricity growth. 
Finally, binary stars typically have much larger mass ratios. Once q > 0.2, the 
3:1 exciting resonance no longer lies within the stable disc, and eccentricity 
growth is suppressed. Although work on superhumps invariably uses a smaller 
mass ratio (to place the resonance within the disc), other resonances from the 
secondary will be far stronger in such calculations, and it is not clear what 
their effect will be. 



The growth of eccentricity demonstrated in this paper has important impli- 
cations for the possibility of hole growth in protoplanetary discs. Accretion 
discs containing holes have been observed in t hree systems: CoKu Tau/4 



(iD'Alessio et all l2005h . GM Aur and DM T au dCalvet et all . I2005T ). It has 



already been suggested (jOuillen et all 12004 ) that CoKu Tau/4 contains a 
planet. Other models have been sug gested, for example the photoevaporation 
model of Alexander et al.l (l2006al ibl) . and there are presently insufficient data 
to rule all of these out. However, in CoKu Tau/4 at least, the edge of the 
accretion disc is observed to be very sharp (from SED fitting), which lends 
itself to interpretation as the edge of a gap induced by a planet. The 1000- 
2000 orbit timescale found in the calculation of section [2] fits easily within the 
estimated ages of CoKu Tau/4, GM Aur and DM Tau. 



Consider the case of CoKu Tau/4: CoKu Tau/4 is a comparatively young 
(~ 1 Myr), low mass (M* ~ 0.5 M Q ) system with a large inner hole in its 
accretion disc. This hole stretches out to a radius of ~ 10 au, and is depleted 
by a factor of about 10 5 in dust. There is no accretion signature, suggesting 
that little gas is present in the hole. Numerical e xperiments find that oute r 
disc lies a couple of Hill radii from the planet (e.g. Ide Val-Borro et al.l . 120061 ) . 
so a planet with a mass ratio 10~ 3 would be in orbit at around 8 au, with an 
orbital timescale of approximately 30 yr. Two thousand orbits would pass in 
less than 10 5 yr, comfortably less than the age of the system. However, the 
formation time of the planet, plus the extra time necessary to produce a 10~ 5 
depletion in surface density (around ten exponential depletion times) should 
be borne in mind. 



The constraints are less severe for GM Aur. This is a 1.2 M Q star, which is 
still accreting some gas. There is an optically thick outer disc, truncated at 
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24 au. However, there is extra emission from the inner regions, co nsistent with 



an op tically thin inner dust disc, extending out to around 5 au. ISimon et al. 



( 120001 ) estimate the age as 3 Myr. If the disc edge is being maintained by a 
10~ 3 mass ratio planet, we would expect it to be located at about 22 au, with 
an orbital period of 94 yr. The presence of the optically thin inner disc would 
be consistent with our timescale. If we assume that the planet will have formed 
in around 1 Myr, ten such timescales remain with the age of the system. The 
eccentricity gro wth rat e is exp ected to scale with the square of the mass ratio 
(equation 66 of iLubowl . Il991al ) . Putting this into the timescale constraints, we 
conclude planet with only one third the mass ratio shown above (i.e. 5 x 10 -4 ) 
would probably be sufficient to produce this system. 

It is possible that we are observing an evolutionary sequence, with GM Aur 
being the 'youngest' system and CoKu Tau/4 being the 'oldest'. We place 
'oldest' and 'youngest' in quotes, since they refer to number of orbits completed 
by the planet, rather than the absolute age. Further observati ons of these 



system s will be useful in constraining their properties. Recently, Eisner et al. 



(120061 ) managed to resolve the inner edge of the accretion disc hole in the TW 
Hydrae system. Hopefully, similar techniques could be used on these three 
young systems. 



5 Conclusion 



In this paper, we have reported a new mechanism for producing a hole in an 
accretion disc: eccentricity excitation by a protoplanet. Scaling our results, 
we have shown that the holes in the discs of CoKu Tau/4, DM Aur and GM 
Tau can all be explained by the presence of a ~ 1 Mj planet in these systems. 
Future work should analyse the clearing mechanism in more detail, exploring 
the effects of varying viscosity and planet masses. Performing a parallel study 
with a separate code will be extremely useful. Although we have no reason 



to bel ieve that Fargo is not performing well, the work of Ide Val-Borro et al 



(120061 ) underlines how different codes can give subtly different answers to the 



'same' problem. 

In gener al, we expect a hole to form so long as a gap can form (q > A0TZ~ 



see e.g. iBryden et all Il999l ) and the disc is insufficiently massive to cause 
significant migration (roughly, M disc < M p ). We also require that the viscous 
timescale is longer than the timescale for eccentricity growth (this requires a < 
10~ 2 ), or viscous clearing will dominate. If the gap opening criterion is fulfilled, 
the planet will then prevent the outer disc replenishing the inner, regardless of 
the mechanism which clears the inner disc. In this paper, we have concentrated 
on the growth of eccentricity, but viscous clearing and Lindblad waves are 
alternatives. We shall explore these in future work. The recent improvements 
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in observations provided by Spitzer open up exciting possibilities in the field 
of planet formation. There is now a real possibility that we will be able to 
observe planetary systems in the process of formation. This will provide strong 
constraints on theoretical and numerical predictions such as ours. 
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